CAPTAIN visualisation workflow

Author

Théophile L. Mouton

Published

September 24, 2024

Visualise CAPTAIN results

R libraries

Code 
library(readr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(dplyr)
library(gridExtra)
library(biscale)
library(colorspace)
library(grid)

Budget 0.1

EDGE2

Code 
# Buget 0.1
# Read the RDS file from the Data folder
data <- readRDS(here::here("Data/EDGE_full_results_averaged_budget0.1_replicates100.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the data to sf object and project
data_sf <- st_as_sf(data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create the outline from the projected world map
outline <- st_union(world_projected) %>% st_boundary()

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create the plot
ggplot() +
  # Add points for each planning unit, colored by priority
  geom_sf(data = data_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  # Add the world map with a light gray fill
  geom_sf(data = world_projected, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  # Add the black outline around the globe
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +  # Add the globe border
  # Use a white to yellow to blue color gradient
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                           title.position = "top", title.hjust = 0.5)
  ) +
  # Add labels and title
  labs(title = "Global Distribution of Conservation Priorities",
       subtitle = "Index: EDGE2, Budget: 0.1, Replicates: 100",
       x = "Longitude", y = "Latitude") +
  # Adjust theme
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

Code 
# Save the plot
#ggsave("priority_world_map_mcbryde_outline_EDGE2_0.1_100reps.png", width = 15, height = 10, dpi = 300)

# Protection fraction summary
# Load required libraries
library(ggplot2)
library(jsonlite)
library(here)
library(gridExtra)

# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE_0.1.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))

# Extract Species and EDGE2 from sp
Species <- sp$EDGE2$info$Species
EDGE2 <- sp$EDGE2$info$EDGE2

# Add Species and EDGE2 to prot_frac
prot_frac$Species <- Species
prot_frac$EDGE2 <- as.numeric(EDGE2)

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for EDGE2
hist_edge2 <- ggplot(prot_frac, aes(x = EDGE2)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of EDGE2 Scores",
       x = "EDGE2 Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = EDGE2, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  labs(title = "Scatterplot: EDGE2 vs Mean Protect Fraction",
       x = "EDGE2 Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_edge2, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code 
# Optionally, you can save the grid plot
# ggsave("grid_plot.png", grid_plot, width = 12, height = 10)

FUSE

Code 
# Read the RDS file from the Data folder
data <- readRDS(here::here("Data/FUSE_full_results_averaged_budget0.1_replicates100.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the data to sf object and project
data_sf <- st_as_sf(data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create the outline from the projected world map
outline <- st_union(world_projected) %>% st_boundary()

# Create the plot
ggplot() +
  # Add points for each planning unit, colored by priority
  geom_sf(data = data_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  # Add the world map with a light gray fill
  geom_sf(data = world_projected, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  # Add the black outline around the globe
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +  # Add the globe border
  # Use a white to yellow to blue color gradient
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                           title.position = "top", title.hjust = 0.5)
  ) +
  # Add labels and title
  labs(title = "Global Distribution of Conservation Priorities",
       subtitle = "Index: FUSE, Budget: 0.1, Replicates: 100",
       x = NULL, y = NULL) +
  # Adjust theme
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

Code 
# Protection fraction summary
# Load required libraries
library(ggplot2)
library(jsonlite)
library(here)
library(gridExtra)

# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_01.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))

# Extract Species and EDGE2 from sp
Species <- sp$FUSE$info$Species
FUSE <- sp$FUSE$info$FUSE

# Add Species and EDGE2 to prot_frac
prot_frac$Species <- Species
prot_frac$FUSE <- as.numeric(FUSE)

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for EDGE2
hist_fuse <- ggplot(prot_frac, aes(x = FUSE)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of FUSE Scores",
       x = "FUSE Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = FUSE, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  labs(title = "Scatterplot: FUSE vs Mean Protect Fraction",
       x = "FUSE Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_fuse, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

EDGE2 and FUSE

Code 
# RGB map ----
# Read the RDS files from the Data folder
edge_data <- readRDS(here::here("Data/EDGE_full_results_averaged_budget0.1_replicates100.rds"))
fuse_data <- readRDS(here::here("Data/FUSE_full_results_averaged_budget0.1_replicates100.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Combine EDGE and FUSE data
combined_data <- edge_data %>%
  rename(EDGE_Priority = Priority) %>%
  left_join(fuse_data %>% rename(FUSE_Priority = Priority),
            by = c("Longitude", "Latitude"))

# Normalize priorities to 0-1 range
combined_data <- combined_data %>%
  mutate(
    EDGE_Priority_Norm = (EDGE_Priority - min(EDGE_Priority)) / (max(EDGE_Priority) - min(EDGE_Priority)),
    FUSE_Priority_Norm = (FUSE_Priority - min(FUSE_Priority)) / (max(FUSE_Priority) - min(FUSE_Priority))
  )

# Transform the data to sf object and project
data_sf <- st_as_sf(combined_data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create the color palette
map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)

map_pal_mtx <- matrix(map_pal_raw, nrow = 4, ncol = 4)
map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
map_pal_mtx[1, 1] <- '#ffffee'

map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))

# Create a custom color mapping function
get_color <- function(edge, fuse) {
  edge_class <- cut(edge, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  fuse_class <- cut(fuse, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  return(map_pal[(as.numeric(fuse_class)-1)*4 + as.numeric(edge_class)])
}

# Apply the function to get colors
data_sf$new_color <- mapply(get_color, data_sf$EDGE_Priority_Norm, data_sf$FUSE_Priority_Norm)

# Create the main plot
main_plot <- ggplot() +
  geom_sf(data = data_sf, aes(color = new_color), size = 0.1, alpha = 1) +
  geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
  geom_sf(data = globe_border, fill = NA, color = "grey70", size = 0.5) +  # Add the globe border
  scale_color_identity() +
  coord_sf() +
  theme_minimal() +
  labs(title = "Bivariate Map of EDGE and FUSE Priorities",
       x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(hjust = 0.5))

# Create the legend
legend_plot <- bi_legend(pal = map_pal, dim = 4,
                         xlab = 'EDGE2',
                         ylab = 'FUSE')

# Define the layout matrix
layout <- rbind(c(1, 1),
                c(2, 3))

# Combine the main plot and legend using grid.arrange
combined_plot <- grid.arrange(
  main_plot, 
  legend_plot,
  rectGrob(gp = gpar(col = "white")),  # Blank space
  layout_matrix = layout,
  widths = c(0.25, 0.75),  # Legend takes 15% width, main plot and blank space take 85%
  heights = c(0.8, 0.2)  # Main plot takes 80% of height, legend row takes 20%
)

Code 
# Display the combined plot
#print(combined_plot)

# Save the combined plot
#ggsave("bivariate_priority_map_with_legend_100reps_budget01.png", combined_plot, width = 15, height = 12, dpi = 300)

Budget 0.3

EDGE2

Code 
# Buget 0.3
# Read the RDS file from the Data folder
data <- readRDS(here::here("Data/EDGE_full_results_averaged_budget0.3_replicates30.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the data to sf object and project
data_sf <- st_as_sf(data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create the outline from the projected world map
outline <- st_union(world_projected) %>% st_boundary()

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create the plot
ggplot() +
  # Add points for each planning unit, colored by priority
  geom_sf(data = data_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  # Add the world map with a light gray fill
  geom_sf(data = world_projected, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  # Add the black outline around the globe
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +  # Add the globe border
  # Use a white to yellow to blue color gradient
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                           title.position = "top", title.hjust = 0.5)
  ) +
  # Add labels and title
  labs(title = "Global Distribution of Conservation Priorities",
       subtitle = "Index: EDGE2, Budget: 0.3 Replicates: 30",
       x = "Longitude", y = "Latitude") +
  # Adjust theme
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

Code 
# Save the plot
#ggsave("priority_world_map_mcbryde_outline_EDGE2_0.3_30reps.png", width = 15, height = 10, dpi = 300)

# Protection fraction summary
# Load required libraries
library(ggplot2)
library(jsonlite)
library(here)
library(gridExtra)

# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE_03.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))

# Extract Species and EDGE2 from sp
Species <- sp$EDGE2$info$Species
EDGE2 <- sp$EDGE2$info$EDGE2

# Add Species and EDGE2 to prot_frac
prot_frac$Species <- Species
prot_frac$EDGE2 <- as.numeric(EDGE2)

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for EDGE2
hist_edge2 <- ggplot(prot_frac, aes(x = EDGE2)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of EDGE2 Scores",
       x = "EDGE2 Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = EDGE2, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  labs(title = "Scatterplot: EDGE2 vs Mean Protect Fraction",
       x = "EDGE2 Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_edge2, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

FUSE

Code 
# Buget 0.3
# Read the RDS file from the Data folder
data <- readRDS(here::here("Data/FUSE_full_results_averaged_budget0.3_replicates56.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the data to sf object and project
data_sf <- st_as_sf(data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create the outline from the projected world map
outline <- st_union(world_projected) %>% st_boundary()

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create the plot
ggplot() +
  # Add points for each planning unit, colored by priority
  geom_sf(data = data_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  # Add the world map with a light gray fill
  geom_sf(data = world_projected, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  # Add the black outline around the globe
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +  # Add the globe border
  # Use a white to yellow to blue color gradient
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                           title.position = "top", title.hjust = 0.5)
  ) +
  # Add labels and title
  labs(title = "Global Distribution of Conservation Priorities",
       subtitle = "Index: FUSE, Budget: 0.3 Replicates: 30",
       x = "Longitude", y = "Latitude") +
  # Adjust theme
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

Code 
# Save the plot
#ggsave("priority_world_map_mcbryde_outline_FUSE_0.3_30reps.png", width = 15, height = 10, dpi = 300)

# Protection fraction summary
# Load required libraries
library(ggplot2)
library(jsonlite)
library(here)
library(gridExtra)

# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE_03.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))

# Extract Species and FUSE from sp
Species <- sp$FUSE$info$Species
FUSE <- sp$FUSE$info$FUSE

# Add Species and FUSE to prot_frac
prot_frac$Species <- Species
prot_frac$FUSE <- as.numeric(FUSE)

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for FUSE
hist_FUSE <- ggplot(prot_frac, aes(x = FUSE)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of FUSE Scores",
       x = "FUSE Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = FUSE, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  labs(title = "Scatterplot: FUSE vs Mean Protect Fraction",
       x = "FUSE Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_FUSE, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

EDGE2 and FUSE

Code 
# RGB map ----
# Read the RDS files from the Data folder
edge_data <- readRDS(here::here("Data/EDGE_full_results_averaged_budget0.3_replicates30.rds"))
fuse_data <- readRDS(here::here("Data/FUSE_full_results_averaged_budget0.3_replicates56.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Combine EDGE and FUSE data
combined_data <- edge_data %>%
  rename(EDGE_Priority = Priority) %>%
  left_join(fuse_data %>% rename(FUSE_Priority = Priority),
            by = c("Longitude", "Latitude"))

# Normalize priorities to 0-1 range
combined_data <- combined_data %>%
  mutate(
    EDGE_Priority_Norm = (EDGE_Priority - min(EDGE_Priority)) / (max(EDGE_Priority) - min(EDGE_Priority)),
    FUSE_Priority_Norm = (FUSE_Priority - min(FUSE_Priority)) / (max(FUSE_Priority) - min(FUSE_Priority))
  )

# Transform the data to sf object and project
data_sf <- st_as_sf(combined_data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create the color palette
map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)

map_pal_mtx <- matrix(map_pal_raw, nrow = 4, ncol = 4)
map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
map_pal_mtx[1, 1] <- '#ffffee'

map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))

# Create a custom color mapping function
get_color <- function(edge, fuse) {
  edge_class <- cut(edge, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  fuse_class <- cut(fuse, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  return(map_pal[(as.numeric(fuse_class)-1)*4 + as.numeric(edge_class)])
}

# Apply the function to get colors
data_sf$new_color <- mapply(get_color, data_sf$EDGE_Priority_Norm, data_sf$FUSE_Priority_Norm)

# Create the main plot
main_plot <- ggplot() +
  geom_sf(data = data_sf, aes(color = new_color), size = 0.1, alpha = 1) +
  geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
  geom_sf(data = globe_border, fill = NA, color = "grey70", size = 0.5) +  # Add the globe border
  scale_color_identity() +
  coord_sf() +
  theme_minimal() +
  labs(title = "Bivariate Map of EDGE and FUSE Priorities",
       x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(hjust = 0.5))

# Create the legend
legend_plot <- bi_legend(pal = map_pal, dim = 4,
                         xlab = 'EDGE2',
                         ylab = 'FUSE')

# Define the layout matrix
layout <- rbind(c(1, 1),
                c(2, 3))

# Combine the main plot and legend using grid.arrange
combined_plot <- grid.arrange(
  main_plot, 
  legend_plot,
  rectGrob(gp = gpar(col = "white")),  # Blank space
  layout_matrix = layout,
  widths = c(0.25, 0.75),  # Legend takes 15% width, main plot and blank space take 85%
  heights = c(0.8, 0.2)  # Main plot takes 80% of height, legend row takes 20%
)

Code 
# Display the combined plot
#print(combined_plot)

# Save the combined plot
#ggsave("bivariate_priority_map_with_legend_100reps_budget01.png", combined_plot, width = 15, height = 12, dpi = 300)